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Errata to the NPO Thesis entitled 


A COMPUTER PROGRAM FOR THE ANALYSIS OF TWO-DIMENSIONAL 
HEAT CONDUCTION USING THE FINITE ELEMENT TECHNIQUE 


by Allan Bischof Chaloupka 
June 1969 


In the main program listing: 
page 72 . Change card number 155 to read 


DO 601 N=1, NUMNP 1995 


page 75 Remove card preceding card number 180 that reads 


189 CALL SYMSOL (1) 


page 74 Change statement number on card number 203 to read 


1020 FORMAT (615, D10.3) 203 


In SUBROUTINE FORM: 


page 76 Insert a missing card between card numbers 289 
and 287 to read 


CV = SPHT(MTYPE) * CFUNC(MIYPE,TMEAN) 7 286 


In the MAIN Program listing: 
page 72 Change Card Number 138 to read 


TEMP=TEMPB(N) - TO 138 


Change the + signs to = signs on the 
following Card Numbers: 


144 Surote: i (ee elie 
145 


146 
147 
149 7 
151 
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I, INTRODUCTION 


In this chapter the reasons and advantages of the finite element 
analysis method are given. The purpose of this presentation and the ob- 
jectives of the study are discussed. 

A, REASONS FOR THE FINITE ELEMENT ANALYSIS 

With the ever increasing complexity of problems occuring from engi- 
neering situations in design and analysis, there is a rising need for an 
approximate method of solution for these problems, This technique must 
be flexible enough to encompass a large variety of problems and still give 
accurate answers to the simplest of situations. Since most approximate 
methods of analysis of complex engineering problems involve many tedious 
calculations, it would enhance the usefulness of the technique if it were 
capable of being programed for computer solution. This would relieve 
the engineer of the tedious burden of hand calculated solutions and leave 
him more time to devote to the job of engineering. 

The finite element method of analysis is an approximate method of 
analyzing engineering problems which meet these specifications, It has 
a high degree of flexibility and can be readily applied to computer sol- 
utions to obtain rapid solutions to complex problems. 

Bo OBJECTIVES OF THIS STUDY 

The major objective of this study is a computer program using the 
finite element method of analysis, for the solution of two-dimensional 
unsteady heat conduction problems, In fulfilling this goal there were 
several minor objectives accomplished along the way. The original pro- 
gram was written for the IBM 7094 system. It was converted for operation 


on the IBM/OS 360 computer at the Naval Postgraduate School Computer Center. 


1l 


Involved in this conversion was the changing of control cards and a few 
FORTRAN IV language differences. Also faster random access disks were 
substituted for the magnetic tapes originally used for temporary storage. 
The input and output formats were overhauled to make better use of space 
on the printed output page. More extensive labels were added to the 
printed output in order to ease the use of the program by persons not 
familiar with the finite element technique. 
C, PURPOSE OF THE PRESENTATION 

The purpose of this presentation is to make available a computer 
program using the finite element method of analysis for the solution of 
unsteady heat conduction problems, In doing so an introduction into the 
formulation of the finite element method is given for those not familiar 
with the technique, to give an indication of the power of this type of 
analysis. A description of the components of the program is included to 
enable the user to better understand the logic of the computer solution 
and increase the user's ability to formulate problems to which the pro- 
gram can be applied. Specific instructions on formulation of problems 
and loading of input data is made available in a chapter on user infor- 
mation. Finally, as is customary in works of this type, recommendations 
for future modifications of the program and possible applications are 


included. 


ilk 





II. THE FINITE ELEMENT FORMULATION 


In this chapter the finite element method is formulated. The idea 
of a functional of a function and the variation of a functional are intro- 
duced, These principles are developed into an approximate solution of 
the two-dimensional unsteady heat conduction equation using a matrix 


iterative solution. 


A, FORMULATION OF THE VARIATIONAL METHOD 
If a functional of a temperature field can be found such that its 


extremum yields an equation describing unsteady heat conduction, 


Ri Ty = ‘en (2-1) 
and an equation describing the heat flux across the surface of the region, 
ni key f me ei (2-2) 


then the variation of the functional will yield a set of first order, 
ordinary differential equations in terms of the nodal temperatures, Using 
the method of tensor notation, a comma denotes differentiation with res- 
pect to the following subscript and repeated subscripts imply summation. 
Let JT be a functional of the temperature field T(x,y,z,t) and the 


first time derivative of the temperature field T(x,y,2,t), be defined by 
TCT, T)= (zh Ry reel T ZI )dv + | n.aTas (2-3) 


' e -) 
The variation of TI(T,T) with respect to T (with T held constant) 


is given by 
dt= Si(r+ éA,T) (2-4) 
Se 
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where € is a small parameter and A is one of a family of functions. 
A is zero on the boundary of the region and arbitrary elsewhere. An 


extremum of the functional TI(T,T) implies that SIT, T) equals zero, i.e., 


Str (T+€A,T) =A (2-5) 
JE €=0 


First 


TUTHEATI= | (E(THea): Ris(TEAD 
Vv 


+ec(T+€A)T- E(T +62) Jay nega A) dS (2-6) 
S 


Using the identity 
oy Led pay T+ 6) j= tr. Or)}ae Ry ~ Res (Tt € rit r 
ge 7 


if the conductivity tensor is symmetric i.e. kijj = kji. Then 


6m (r+éA,T) = \(Er-eac asa a. Ry; (7 EN, A 
de y | = 


tecrT -Gd)dv - | nige Aas (2-7) 


The volume incegrat (Ty Re5 A) sd can be transformed into a surface 
V 


integral, 


\ (1; Rij ),s dv = ni kis 1,7 Ads (2-8) 


y S 
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When (2-7) is evaluated at EG =0 


suct +E T) =\Ge Rej Ns; = Reg yeh 


+ec\T Ea) ay -\ nide Ads 
3 


dt(T+€A,T) E CRa Th: v ect-4) Ady 
tae Wel Rey yj Ads - AUGi Ads (2-9) 
$ = 
The extremum of IT ( St 0 )reduces to 


Rg Tos; = #8) —G (2-10) 


in the region V. Also on the surface § 


neRij 1,5 = he Vi (2-11) 


L 


These two equations (2-10, 2-11) are the unsteady heat conduction equation 
and the boundary flux equation. This verifies that a function T which 
yields an extremum of the functional defined by (2-3) satisfies both the 
unsteady heat conduction equation in the region V and the boundary flux 


equation on the surface $ [e}*. 


B, FORMULATION OF THE FINITE ELEMENT PROBLEM 

If the choice of T and T is restricted so that they are represented 
by certain constants which are obtained in their formulation, the func- 
tional [{f becomes a real-valued function. In this case TI shall be 
Tr} »{TT) where {tT} and (T} are vectors of nodal point values 


of {r} and {72 The extremum of T(r > 7} ) must now be found which 





*Numbers in brackets refer to similarly numbered items in Bibli- 
ography. 
15 





implies 


Sir 72, (7) =) (2~12) 
OT; 


From this point on the finite element formulation will be considered 
in two dimensions only. The region V will be divided into a number plane 
triangular elements. The functions T and T can be uniquely described 
within each of these elements by the values of the function at each of 
the nodal points of the element Tj, Tj, Tm, and linear function of the 
coordinates of the element. 


Let the temperature field in an element be given by 


TL, 4) = Boy fAGO, F=43,m aay 


and the time rate of temperature change given by 


TY, 0) = GO, W>[AlL To a 


<5> is a vector which specifies the spacial approximation of T and . 
The matrix 4] is a constant and is defined by the above relationships. 
Its exact value and derivation will be discussed in detail in section II, 
ere 

The temperature gradient Tan Cn “ x,y) can be expressed in terms of 


the nodal temperatures by 


Le = 
Ty DO wity =a a 


In the following development the subscript ; , denoting nodal point 


values, will be dropped and will be assumed. 


16 


Writing the functional TToin terms of nodal “= values, 


TTL (EMAC TIPE ciples oan 


— {TY AT @y'G, dy - | Tale, Was — 
Taking the variation of VP with respect to iT} and setting it equal to 

an itt $4}) = IX {Tt | -{e4} }+{c] 7 34 =0 (2-17) 
where 


li 


LK] 
Ic] 


| ay (oy IRI IfAl dv (2-18) 


{i 


ec lal) > [Al av _ 
|e" = \ % Al @> g Vv (2-20) 


er} = | BT @Find feos sea 


C. BOUNDARY CONDITIONS 

There are four types of boundary conditions to be considered. (1) 
specified temperature, (2) specified flux, (3) convection boundary layer, 
and (4) adiabatic, For a specified temperature boundary condition T,; is 
constant over boundary segment 4 For a specified flux q = q over 
boundary segment S9, For a convection boundary layer q = h(Ty - To) over 
boundary segment 53. Finally for an adiabatic boundary condition q = 0 


over boundary segment 54, 
: 17 


The boundary integral (2-21) becomes 


a“, 

oe 

at 
ji 


ag af \ui¢Tt =e (2-22) 
| [Al int {L$ ds (2-23) 


i 


where {ox} 


TH) =k (A) oA as (2-24 


“AG 


) | i 
yh he (Ale) ds (2-25) 
Ss 
The integral on segment Sy is zero because the variation of the functional 
was specified as zero on the boundary. The integral on segment Sy, is 
zero since there is no heat flow across the boundary. 
To obtain the extremum of the functional {[ the nodal point tempera-~ 


tures (Tt which satisfy the following matrix equation must be found [6]. 


((KI- HI rt + Kedtas - \% *| ~\ 4 + Shri = © (9296) 


D, THE MATRIX EQUATION 
Let the time variable be tj such that tj = iat i ="O7), 2a. 


For this development i will refer to time increments. 


\K| ~ LK |-[H] (2-27) 
Equation (2-17) can be written as 


[K] 1} a Ic] (TS in {ft (2-28) 
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Let 


where = {3" te } - -(h"} ide 


e6 ; 
Assuming that . is constant over the interval ty < tet, +p Lm@., 


Th, = ({TI a - {Th Yat 


a). ; can be obtained by a Taylor's expansion around t = ty . 
«{? 


re, oh wails stk — 
(Tha, ATE. atlThe + SETA) aan 
Wi. = S7i, + 2 (ET... 7 {T2.) (2-32) 


Using equations (2-29) and (2-32) we and iM can be determined 


interms of {T{, and {Tf}, . Solving (2-32) for {Tf}, gives 
rhe, 7 a Thi, —{T}.) — {Th (2-33) 
Substituting this into (2-28) gives 
rR. Fhe + AGT 7h) aw 
Also from equation (2-28) 


rt. = fel f}- ©] is: "= 


Substituting (2-34) into (2-35) gives 
(ii le +5 2 [a)) The, =(Alel kr}. + _ + ft (2-36) 


((K] (4) (AL = we (OTT h(t als Buf.) ean 


In a multi-element body the element matrix equations must be assembled 
into a single matrix equation. The assembly and solution of the equation 


is a complex task but it can be greatly simplified by the use of a com- 


puter [6]. 


E, FORMULATION OF THE ELEMENT MATRIX 
If the region to be analyzed is divided into plane triangular elements 


then a typical element would appear like the following figure. 





Figure 2-1 Triangular Element 


The temperature in an element can be described in terms of a linear com- 
bination of the coordinates of the elements nodal points and the values 
of the temperature at each of the nodal points. 


Assuming that the temperature distribution is of the form 


TOGY,£) = & W4a, ax +4, (2-38) 
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The temperature at each of the nodal points can be written as 


Ti = A, + XA “a Ye Ag 
Vy = A + | i } ¥y X 3 (2-39) 


en - ww + X vate ¥ Yn A3 
Solving for the coefficients Ai Aas and Ax, using Cramer's rule the 


following results are obtained. 


1 Ke VF 
Ler Xy 1, = £A (2-40) 


Note: is the area of the triangular element. 


Ti XY L 
A=|T x ¥% |/ 44 
Tw Amn Yn 


‘mete, 
—} 


ap 
Sa LA 
lb Tae eR 
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x 
A; 
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yvi 


—_ 
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el 
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Zi, 


aaa i [te Ym | + lees 2 ¥{) 2h 


Y, \ vi 
“+ To 1) ZA 


My] _ dt Xm 
“le My red il 


me em) 4k 
= (TPH 
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l, 
mI 


























Tee | 
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Recalling that the temperature distribution was previously described by 


the expression 
= <P>|A\{T} (2-42) 


+7 is the vector of the nodal temperatures and \A] is the spacial 


dependence matrix which is a constant. let 


<d> =3 a BXBY> (2-43) 


a Vy A12 Hz 
‘Th = 4% IAL = Jay a; 
eli A 3 A352 Ase 


then 


T = <P>[alor} 


Ze 


| Gu li+ of, | + oe 
| t ~ </ x y } (4, T; + a], + es low) (2-43a) 
\.. K a X55 [i+ a5 hn) 


T= (48.7; +03 Ta)+ Cal thal HOt) x +@%,Te 4a 52 1) 4Vsg lin) ¥ 


Therefore, comparing (2-43a) with (2-41) 
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The spacial dependence matrix is then defined by 


nee Xray’) ae a Yate Y fe Ky ye 


[A] = Ly oe 1 eae (2-44) 
ZA 


F, THE METHOD OF LUMPING ELEMENT PROPERTIES 
When working with a computer it is usually more convenient to work 
with a quadrilateral element. The quadrilateral is divided into four 


triangles with the center nodal point common to each of the triangles. 


(see figure 2-2). 





Figure 2-2 Quadrilateral Element 
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The matrix equation (2-28) can be assembled for this element in the 


following form 





Kii Ki2 Ki3 Kiq4 Fis Ci3 C14 S15 

Ko1 Koo Ko3 Koq Kas Co3 C24 C25 

sa Pyarha3 Kay M5 a5 “Gy, “35 
41 ®42 *a3 Sag Kas C43 [44 [45 
Ba oes gh Noy, #55 | C5g0'st “s5q 


Partitioning this equation results in 
Kei] {Kis} {r,} (a5) §% 7 ty f. 
i i5 i i5 
[3 yh + ; i = fe (2-45) 
Ks i> ¥55 i) (653) Css be ees 


i,§e: 1,272,394 


This equation can be multiplied out to yield two equations 


Ki \{ts + {Kist ae [¢.5] Sty + 4C15¢ i. XS J £34 (2-46) 
(Ks) tah + Bom tf CAPS Ess + Gaels ter Ge (247) 


The specific heat matrix [c ] is approximated by lumping the heat ca- 
pacities at the four external nodes. Thus jc] becomes a diagonal matrix 


(4x4) and C55 = 0. Therefore (2-46) can be written as 


[Ki 5] {7,2 + Viera To [ca] t,t " {st (2-48) 
and (2-47) as 
(55) 17,4 + Bee's = Us (2-49) 


Solving (2-49) for T, and substituting into (2-48) results in 
= 
[Ki] Stuf + {Kis} Kgs (f5 - (Ks 5)(T:\ ) + [cis] { - 5 £4 (2-50) 


25 


( (kK. J = Kt Ke cK.) ) ir} + fess} rs ff} al {Kis} Kea (2251) 


Equation (2-51) is analogous to equation (2-28) except that [K] and Ic] 


are now (4x4) matrices and tf I is a (4x4) vector (6 } Z 
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III, PROGRAM STRUCTURE 


This chapter describes the Aeneaail composition of the program and 
gives a description of the function of each section of the program. 
A, PROGRAM IDENTIFICATION 

The program was written by Robert E, Nickell while in the Department of 
Civil Engineering, University of California at Berkeley. It was formerly 
identified as W2498,W035--2--57343, Nickell. It was written in July 1968 
and was revised by Allan B. Chaloupka in May 1969. Formulation was 
based upon a variational principle derived by M. E. Gurtin of Brown 
University (4). 

The program has been renamed DFETHC-Finite Element Transient Heat 
Conduction. 
B, PURPOSE 

The purpose of the program is to produce high speed solutions to 
transient heat conduction problems, These problems may involve plane or 
axisymmetric geometry, non-~homogeneous anisotropic material configura- 
tions, temperature dependent material properties, and time dependent 
boundary conditions. 
C, PROGRAMING INFORMATION 

The program was originally written in FORTRAN IV for the IBM 7094 
computer, It was converted for use on the IBM 0S/360 Model 67 computer. 
The program was also rewritten in double precision (REAL *8) , 
D. PROGRAM CAPACITY 

The size problem that can be analyzed by this program is subject to 
the following limitations" 


Massiaam nmumberWef medal pedmts . . cigs « cles ile cic ccc tive ecce JOO 


Zu 


Maximum number of quadrilateral elements ......c.ccccceccececcse 490 
Maximum number Of Mate mila lism acs «ss 66 0050 000 os 6 opclls cercte sn 
Maximum difference of nodal point numbers 

for the samerelemenmtr a, .... sss se +6 66.0. ote oicie.c.0 ¢ictteneeets/ sane 

These are not stringent limitations and they may be adjusted some- 
what by the user to fit individual problems. The limits are changed by 
changing the appropriate DIMENSION and COM’: statements and one card in 
the segment that calculates the half-band width, The only real limita- 
tion placed upon the size of the problem that can be analyzed is the 
amount of core storage available in the computer used. 

E, GENERAL STRUCTURE 

The program consists of a main program and several subprograms used 
for repetitive operations. The main program can be divided into three 
parts. The first part is concerned with reading into the program and 
printing out the data describing each problem. It is portrayed by the 
flow diagram in figure 3-l. 

The second portion of the main program is concerned with the routing 
of operations according to certain logical variables, These variables 
describe four possible combinations, They are constant material proper- 
ties and constant boundary conditions, temperature dependent material 
properties and constant boundary conditions, constant material proper- 
ties and time dependent boundary conditions, and temperature dependent 
material properties and time dependent boundary conditions. Each of these 
possibilities is portrayed by a flow diagram in figures 3-2 through 3-5. 

The third portion of the main program is concerned with calculating 
the effective forcing vector, solving the matrix equation and printing 


the temperatures for each time increment. The operation is then directed 


ze 
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Figure 3-1 Main Program Flow Diagram, Part I 
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(GAPE SYMSOL(1)  ) 
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Figure 3-2 Main Program Flow Diagram, Part II 
Constant Properties and Boundary Conditions 
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Figure 3-3 Main Program Flow Diagram, Part II 
Temperature Dependent Properties and 
Constant Boundary Conditions 
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Figure 3-4 Main Program Flow Diagram, Part Il 
Constant Properties and Time Dependent Boundary Conditions 
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Figure 3-5 Main Program Flow Diagram, Part II 
Temperature Dependent Properties and 
Time Dependent Boundary Conditions 
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back to the beginning of the second section of the main program to start 
calculations for the next time increment, The third portion of the main 
program is portrayed by a flow diagram in figure 3-6. 
F, MAIN PROGRAM 

The main program is composed of a number of individual steps which 
comprise the finite element solution technique, and a system for the in- 
put and output of problem information. Not included in the main program 
are the matrix assembly processes, the matrix equation solving processes, 
and the boundary condition input routines, These are included in sub- 
programs FORM, SYMSOL, TFUNC, CFUNC,FCBC,FIBC, and FFBC, 

1. Control Information 

The first information to be read into the program is the control 
information, This information sets the basic variables for each of the 
steps of the program and also routes the problem to various sections of 
the program. This information includes the number of nodal points, the 
number of elements, the number of convection boundary conditions, the num- 
ber of materials, the number of time sequences, whether the problem is in 
plane or axisymmetric geometry, the output print interval, the system of 
units, whether a final spacial temperature distribution is desired, and 
the reference temperature of the body. 

Also read in with the control information are two titles, The 
first is a 72 space title that can be used by the user to label each 
problem. The second is a 48 space title that is used to input the labels 
used with the system of units chosen. 

The appropriate titles and control information are printed on the 


first page of the printed output, for each problem, 
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Figure 3-6 Main Program Flow Diagram, Part III 
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2. Material Properties 


The thermal conductivity, density, specific heat, and heat gener- 
ation properties for each material are read into the program, The materi- 
al properties are then printed in order on the printed output page and 
labeled with the appropriate units. Since the program is for anisotropic 
materials, Kxx»Kyys and Kyy are read in order to form the thermal con- 
ductivity tensor. The specific heat and density are read in as a product 
in order to reduce the number of variables in the program. 

3. Time Sequencing Information 

In this section the time sequencing information is read in for 
each sequencing scheme, This information consists of two logical varia- 
bles and two numerical variables. The logical variables indicate whether 
there are temperature dependent properties and if there are time dependent 
boundary conditions. The two remaining variables describe the number 
of time steps and the size of each time increment. The time sequencing 
information is printed on the printed output page with appropriate labels. 

4. Nodal Point Information 

Ordinarily this segment of the program will read the x,y coordinates 
and initial temperature of each nodal point and print the information 
with the appropriate dimensions. Because of the possibility of a large 
number of nodal points (500) and the equally large possibility of human 
error in punching data cards there is a segment of the program that will 
generate nodal point information. This segment has a few limitations, 
which will be discussed in chapter IV, but it is usually able to generate 
almost all of the interior nodal points. As an example in the configura- 
tion shown in figure 3-7 the nodal points indicated by squares were read 
into the program and the remaining nodal points were generated by the 


program. 
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Figure 3-/ 





Nodal Point and Element Generation 
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5. Element Information 

The purpose of this segment is about the same as the previous 
one, The number of the element, the number of each of its nodes, the 
type of material, and the amount of heat generation within the element 
are read into the program and printed out with appropriate labels. Again 
this segment has the capability to facilitate the loading of data by 
generating element information. This segment also has some limitations 
which will be discussed in chapter IV, but it is usually able to generate 
most of the elements in a given configuration. As an example, in figure 
3-7 all the elements were generated except those that are starred. 

Also included in this segment is a routine that calculates the 
half-band width of the problem. This is an important control parameter 
that is used later in the program to assemble the matrix equation. 

6. Time Steps 

The time step portion of the program includes two large loops. 
The first loop executes the program through each time sequence of the 
problem. The second loop executes the program through the required num- 
ber of time increments for each time sequence. Inside these two loops 
the two logical time sequencing variables and the loop counters control 
the routing of the problen. 

The two logical variables control the choice of variable material 
properties and boundary conditions. If the material properties and bound- 
ary conditions are constant then the coefficient matrix and forcing vector 
are calculated only once, The matrix equation is then solved for suc- 
cessive time increments using only the original matrices. If variable 
material properties are involved the effective coefficient matrix must 
be reapplied since the coefficient matrix is reformed from zero values. 


If time dependent boundary conditions are involved then the contribution 
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of the boundary conditions to the coefficient matrix must be reapplied. 
The coefficient matrix from the original material properties is kept in 
temporary storage. When both time dependent boundary conditions and 
temperature dependent properties are used then the coefficient matrix 

and forcing vectors must be reformed for each time increment based upon 
the immediate time and temperature values and initial properties. The 
data for temperature dependent properties is calculated from the initial 
conditions and the present nodal temperatures, The data for time depend- 
ent boundary conditions is either read into the program with each time 
increment or generated within the program. 

Different time sequences can be run for the same problem, This 
enables the user to speed or slow the time change during a particular 
segment of the problem. 

7. Convection Boundary Conditions 

A convection boundary condition occurs when the problem includes 
a convection boundary layer as part of its boundary conditions. One con- 
vection boundary condition consists of two nodal point numbers, the con- 
vective heat transfer coefficient and the ambient fluid temperature, The 
two nodal points define the convection boundary for each condition. 

Time varying convection boundary conditions can be obtained by 
varying the ambient fluid temperature, The new fluid temperature may be 
read with each new time step or it may be generated by using the function 
subprogram FCBC, FCBC is used in the same manner as the function sub- 
programs for temperature dependent material properties, The argument of 
FCBC is TIME and the value of FCBC is currently set at one. If tempera- 
ture dependent material properties are used then the convection boundary 
condition section stores and reapplies the convection boundary condition 


contributions to the coefficient matrix and the forcing vector. 
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8, Temperature and Flux Boundary Conditions 


A temperature or flux boundary condition occurs when the problem 
specifies a temperature or flux at a nodal point. One temperature or flux 
boundary conditions consists of the nodal point number and the value of 
the temperature or flux at that nodal point. 

Time varying temperature and flux boundary conditions can be ob- 
tained by varying the temperature or flux at each nodal point, The new 
temperature or flux may be read with each new time step or it may be 
generated by using the function subprograms FIBC or FFBC for time depend- 
ent temperature or flux respectively. FIBC and FFBC are used in the 
Same manner as the function subprograms for temperature dependent materi- 
al properties. The argument of each is TIME, Their current values are 
both one, If temperature dependent properties are used the temperature 
and flux boundary condition section stores and reapplies the boundary 
conditions to the coefficient matrix and the forcing vector. 

9. Temperature Output 

This segment of the program is concerned with the print out of 
the temperature solutions for each nodal point at each time step. The 
value of the time is printed and the nodal point numbers and correspond- 
ing temperatures for that time are printed in rows below. If it is de- 
sired, the temperature can be printed with spacial data at the end of 


each time sequence. 


G. SUBROUTINE FORM 

The subroutine FORM is used to calculate the effective element ther- 
mal conductivity, heat capacity, and heat generation. lhe subroutine 
performs these operations for each of the elements of the configuration 


and, if temperature dependent properties are used, for each time increment. 
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Figure 3-9 Computer Storage for the 
Matrix of Coefficients 
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First the mean x,y coordinates and the mean temperature of the element 
are calculated by averaging the values of each of the nodes. Using the 
mean x,y coordinates the quadrilateral element is divided into four tri- 
angles. The conductivity tensor, heat capacity, and heat generation 
vectors are formed for each of the four triangular sections. As these 
matrices are formed they are assembled into the larger (5) nodal point 
equation for the quadrilateral element. The heat capacity and heat gen- 
eration vectors are assembled so that the fifth term of each vector is 
dropped when forming the lumped properties matrix equation. The center 
nodal point is eliminated from the thermal conductivity tensor (5x5), 
according to equation (2-52), to give the thermal conductivity tensor 
(4x4) for the lumped properties matrix equation. The thermal conductivity 
tensor and heat capacity and heat generation vectors for each element 
are returned to the main program where they are loaded into the coeffi- 
cient matrix and vectors used in the solution of the total matrix equa- 
tion for each time step. For problems where the properties of the 
materials are not temperature dependent the thermal conductivity tensor 
and the heat capacity and heat generation vectors are calculated only at 
the beginning of the problem. If temperature dependent material proper- 
ties are included these matrices are calculated at the beginning of each 


time step. 


H, SUBROUTINE SYMSOL 

The subroutine SYMSOL in conjunction with the loading methods of the 
main program is used to solve the finite element matrix equation. It 
consists of two parts. The first part reduces the coefficient matrix 
to its final form. The second part forms the effective forcing vector 


and solves for the temperatures at each of the nodal points. 
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In an ordinary formulation the final matrix equation for the finite 
element method would look like figure 3-8, The matrix of coefficients 
)R] is a banded symmetric matrix. In the computer, in order to conserve 
core storage, the matrix of coefficients is stored in columns, The first 
column is the diagonal of the matrix of coefficients. (See figure 3-9). 
All the operations performed on the matrix of coefficients are done on 
the upper half of the symmetric banded matrix in its stored form. 

After the coefficient matrix is formed it is put into upper triangular 


form. This yields a matrix equation of the form shown below. 


Ni 


Figure 3-10 Matrix Equation in Upper Triangular Form 


This is accomplished by the first segment of SYMSOL, The second portion 
of SYMSOL, using the effective forcing vector, starts with the temperature 
at the bottom of the nodal temperature vector and back substitutes to 
solve for the nodal temperatures, These operations are performed for each 


time increment. 
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I, TEMPERATURE DEPENDENT PROPERTIES 

There are two function subprograms that are used in conjunction with 
temperature dependent properties. They are TFUNC and CFUNC, TFUNC is 
used for the temperature dependent thermal conductivity and CFUNC is used 
for the temperature dependent product of the density and specific heat. 
In the programs present configuration each of these function subprograms 
is prepared for temperature independent properties and their values are 
therefore equal to one, If the thermal conductivity were temperature 
dependent in the form of: k = k, (1 + bT), where b is a constant coeffi- 
cient and ky is the value of the thermal conductivity at T = 0, then 
TFUNC would have to be modified by the user to include the factor 1 + bT, 


The material properties read into the program would then be zero tempera- 


ture properties. 
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lated. 


IV. USER INFORMATION 


This chapter describes the manner in which problems should be formu- 


included. 


A. 


The first group includes the titles and all the control parameters, 


second group is the material properties. 


INPUT DATA 


quencing information. 


fifth group is element information. 


tions. 


B, 


CONTROL PARAMETERS 


Also the procedures used in inputting data into the program are 


Data for each problem processed by this program is assembled in groups. 


The 
The third group is the time se- 
The 


The fourth group is nodal point information. 


The sixth group is boundary condi- 


After each data deck there must be three blank cards. 


The following control parameters are read into the program. 


NUMNP 
NUMEL 
NUMCBC 
NUMMAT 
NSEQ 
KAT 
INTER 
UNIT 
KPUNCH 
JUMP 
KEY 


TO 


number 


number 


number 


number 


number 


of nodal point 

of elements 

of convection boundary conditions 
of materials 


of time sequences 


type of geometry, O implies plane geometry 


output 


print interval 


units system 


Spacial temperature output 


generated temperature of flux boundary conditions 


generated convection boundary conditions 


initial temperature 
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Also included with the control data are two labels. TITLE (18) is an 
arbitrary label that may be used to identify individual problems. It 1s 
72 spaces long. The second label is UN(12). UN(12) is an array of four 
unit labels. Each has 12 spaces. The units must be printed in the fol- 
lowing order: length, mass, time, and temperature. These labels must be 
left justified. 
©. MeERTAL PROPERTIES 


The following material properties are read into the program. 


XCOND k 
XX 
YCOND Kyy 
XY COND k 
xY 
SPHT ec 
HX heat generation 


Material properties are read in with one material to a card. If the 
maserial is} isotropic k,, = Kyy = 0. If the material is anisotropic the 
three values will not be the same. Each material is given a number ac- 
cording to the order it is read into the program. 

If temperature dependent material properties are used both logical 
variables must be true. Also both JUMP and KEY must be some number other 
than zero and FCBC,FTBC, and FFBC must be one. 

D. TIME SEQUENCING INFORMATION 


The following time sequencing information is read into the program. 


LTAG temperature dependent properties (true or false) 
MTAG time dependent boundary conditions (true or false) 
ITAG number of time steps 

TAG time increment 


The time sequencing information controls the operation.of the program 


through the use of the logical variables LTAG and MTAG but it can also 
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be used to slow or speed time during the problem. An example would be 
two sequences with the first solving for the temperature 60 times during 
the first hour of the problem (i.e., every minute), The second time se- 
quence would solve for the temperature 6 times in the second hour of prob- 
lem time (i.e., every 10 minutes). 
E, NODAL POINT CONSTRUCTION 

The following nodal point information is read into the program for 


nodal point data. 


N number of the nodal point 

KODE l indicates a temperature boundary condition 
XK x coordinate of the nodal point 

ve y coordinate Of the nodal point 

T initial temperature of the nodal point 


The data for each nodal point is put on one card, 

Earlier it was mentioned that the program could generate nodal points. 
This feature has a few limitations that must be observed or incexeean 
information will result. A nodal point whose KODE is 1 cannot be gener- 


ated. Nodal points must be generated along straight lines. The initial 


temperature of the generated nodal points must be the same and equal to IO, 


the initial temperature. To generate nodal points simply leave out the 
data cards between two nodal points on a straight line and the program 
will generate the information for the nodal points in between. 

In constructing the mesh of nodal points for a problem it is sometimes 
impossible to divide the area to be analyzed into uniform rectangular 
quadrilaterals, At times it will not be desirable to do so. Consider 
figure 4-1. In preparing a nodal mesh for this configuration the first 
step would be to assume that by symmetry there is no heat flux across the 


center lines. This reduces the figure to one of its quadrants as in 
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figure 4-2, The size of the quadrilateral elements would vary depending 
upon how interested the user was in a certain area of the figure. In 
figure 4-2 the area around the curve is being more closely examined but 
the mesh size must also be small to enable the user to approximate any 
curved surface using a series of straight line segments. It is almost 
impossible to divide an area with a curved boundary without encountering 
a right triangular element. To convert this to a quadrilateral element 


place a nodal point in the middle of the hypotenuse. 


Figure 4-3 Conversion of a Right-Triangle 
to a Quadrilateral 

In numbering the nodal points it must be kept in mind that the maxi- 
mum difference in nodal point numbers in one element is 31. This is a 
limitation set by the half-band width of the computer program. It is 
convenient to increase the numbering of the nodal points and the elements 
in the same pattern, although it is not a requirement. 
F, ELEMENT CONSTRUCTION 

The following element information is read into the program for ele- 
ment data, 

N number of the element 

IX(N,K) Ke 1-4 numbers of the element's nodal points 

IX(N,5) material type 

HIGEN element heat generation 
The data for each element must be placed on one card. The nodal point 


numbers are placed in clockwise order starting with the number in the 
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upper left hand corner. Elements may be generated by the program but 
there are some limitations. Only elements of the same material type can 
be generated at one time. The heat generation must be the same for the 
elements generated together. Elements must be generated in rows or col- 
umns, To generate elements simply leave out those element cards that 
you wish generated. 
G. BOUNDARY CONDITIONS 

There are two types of boundary conditions that may be read into the 
program, convection and temperature or flux. (An adiabatic boundary con- 
dition is a flux boundary condition with the flux equal to zero.) The 


following information is read into the program for convection boundary 


conditions. 
IB nodal point number 
JB nodal point number 
HB convective heat transfer coefficient 
TEMPB ambient fluid temperature 


The data for each portion of the convection boundary is put one one card. 
If time varying ambient temperature is to be generated by the program the 
control parameter KEY must be some number other than zero, say one. If 
the time varying ambient temperature is not generated by the program new 
data for each segment of the convective boundary must be read into the 
program with each new time step. In this case KEY is zero. 

The following information is read into the program for temperature 
and flux boundary conditions. 

NN nodal point number 

TQ temperature or flux 
The data for the temperature or flux boundary condition at each nodal 


point is put on one card. If time varying temperature or flux is to be 
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generated by the program the control parameter’ JUMP must be some: number 
other than zero. If time varying temperature or flux is not generated 
by the program new data for each boundary nodal point must be read inta 
the program with each new time step. In this case JUMP is zero. Both 
temperature and flux boundary conditions must be introduced into the pro- 
gram in the same manner. For both convection and temperature or flux 
boundary conditions, if time variations are used the logical variable 
MTAG must be true. Also, regardless of the position of the last nodal 
point, a boundary condition card must be read into the program for it. 
If the last nodal point is not a boundary point then the data: card con- 
sists of the nodal point number and a value of zero in the temperature 
or flux field. 
H,.. UNITS 

The units used for each problem must be consistant with the input 
data or the labels on the output will be incorrect. Presently there are 
six different units systems available. Each system consists of a length, 
mass, time, and temperature, Each system has a number which must be read 
in as a control parameter (UNIT), The systems are: 

1. feet, pound mass, hours, °F 

2. inches, pound mass, seconds, °F 

3. feet, pound mass, seconds, °F 

4, meters, kilograms, seconds, °C 

5. centimeters, grams, seconds, °C 

6. feet, pound mass, minutes, OF 

If the user does not prefer the units available he can use another 
system in his input and change the labels or just ignore the labels al- 
together. The system of units will depend greatly upon the configuration 


of the problem and the temperature differences involved. 
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I, ERROR EXITS 

There are three error exits in the program. The first occurs if a 
non-positive number of nodal points is read into the program. This exit 
is used to stop the program, The other two occur in the nodal point and 
element information reading sections. If these statements are executed 
then the execution of the program is terminated, If either of the last 
two exits are executed a statement will be printed to indicate which one 
it was. Usually termination means a bad data card or that the user was 
trying to generate points beyond the capability of the program. 
J, COORDINATES 

The cartesian system of coordinates x and y are used for problems in- 
volving plane geometry. When axisymmetric problems are considered the 
x,0,z system of polar coordinates is used. The analysis of axisymmetric 
problems is conducted in the x,z plane. 
K, TEMPERATURE OUTPUT 

The temperature output can be modified in two ways. The number of 
time increments printed out can be limited by using the control parameter 
INTER, If the program is directed to solve for the temperature every 
minute for 100 minutes and INTER is equal to 25, temperature information 
will only be printed for 25, 50, 75, and 100 minutes. The second manner 
in which the output may be modified is at the end of each time sequence 
the x,y coordinates of each nodal point and the temperature are printed 
if the control parameter KPUNCH equals one. 
L. INPUT DATA PREPARATION 

The following sections describe the manner in which the input data 
is to be placed on data cards, The data is prepared in groups and the 
final data deck is assembled by putting the data groups together in order 


of their listing below. 
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1. Control Information 


(a) Title Card (1844): 72 arbitrary spaces for an alphameric 


problem label punched in columns 1-72 


(b) Units Card (1244): 


Columns 


Labels should be left justified. 


1-12 


13-24 


25-36 


37-48 


Alphameric unit labels 


length 
mass 


time 


temperature 


(c) Control Information Card (1115, F10.0) 


Columns 


1-5 


16-20 
21-25 
26-30 
31-35 
36-40 
41-45 
46-50 
31-55 


56-65 


NUMNP 
NUMEL 


NUMCBC 


NUMMAT 
NSEQ 
KAT 
INTER 
UNIT 
KPUNCH 
JUMP 
KEY 


TO 


number of nodal points. 
number of elements - 
number of convection 
boundary conditions 
number of materials 
number of time sequences 
type of geometry 

print interval 

units system (integer) 
spacial distribution 
temperature conditions 
convection conditions 


initial temperature 


"I' formats must be right justified. 


2. Material Properties 


(a) Material Properties Card (5D10, 3) 


Columns 


1-10 


COND 
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kx 


11-20 YCOND kyy 
21-30 XYCOND Kxyy 
31-40 SPHT ec 
41-50 HX heat generation 
One card for each material. 
3. Time Sequencing Information 
(a) Time Sequencing Card (2L5,110,D10.3): 
Columns 1-5 LT AG temperature dependent 
properties (true or false) 
6-10 MTAG time dependent boundary 
conditions (true or false) 
11-20 ITAG number of time steps 
21-21 TAG time increment 
One card for each time sequence, 
4. Nodal Point Information 


(a) Nodal Point Card (215,3F10.0): 


Columns 1-5 N nodal point number 
6-10 KODE temperature boundary 
condition 
11-20 X x coordinate 
21-30 Bi y coordinate 
31-40 Tk initial temperature 


One card for each nodal point not generated by the program. 
5. Element Information 
(a) Element Card (615,D10. 3): 
Columns 1-5 N element number 


6-10 IX(N,1) nodal point number 
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11-15 IX(N,2) nodal point number 
16-20 IX (N, 3) nodal point number 
21-25 TX(N,4) nodal point number 
26-30 TX(N,5) material number 
31-40 HTGEN element heat generation 
One card for each element not generated. Nodal points are num- 
bered clockwise starting in the upper left hand corner. 


6. Convection Boundary Conditions 


(a) Convection Boundary Card (215,D10.3,F10. 0) 


Columns 1-5 IB nodal point number 
6-10 JB nodal point number 
11-20 HB convective heat transfer 
} coefficient 
, 21-30 TEMPB ambient fluid temperature 


One card for each convection boundary condition. 
| 7. Temperature and Flux Boundary Conditions 
| (a) Temperature or Flux Boundary Card (110,D10.3) 
Columns 1-10 NN nodal point number 
11-20 TQ temperature or flux 
One card for each temperature or flux boundary point. A tempera- 
ture or flux boundary card for the last nodal point must be read even 
if it isn't a boundary point. 
8. Multiple Problems 
The above cards comprise a data deck for one problem. Multiple 
problems may be solved by putting their data decks together, The only 
limitation on this procedure is computer time, It is important that 
there be three blank cards after the last problem data set in the data 


deck. This is to terminate operation of the program. 
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RECOMMENDATIONS 


Since the number of problem types that can be solved by this program 
is limited mostly by the imagination of the user it is highly recommended 
that this engineering tool should be prodern to the attention of those 
students studying in the area of unsteady heat conduction, Only through 
complete familiarization with the program can the user make use of its 
great potential to solve problems involving complex configurations and 
boundary conditions, 

Recommendations for further study in this area are: 

1. Expansion of the finite element method to three dimensional con- 
figurations, 

2. Revision of the boundary condition segments of the program to 
make use of random access disks for temporary storage. 

3. Construction of a supplementary program that will punch input 
data cards for complex time dependent boundary conditions. 

4, An isothermal contour plotting subroutine capable of being used 
in conjunction with the present NPS Computer Center subroutine DRAW, 

5. Use of the program in conjunction with the NPS Computer Center 
IBM Optical Display Unit, in order to observe transient heat conduction 
differences after minor alterations in mesh configuration. 

6, Use of the program in conjunction with existing finite element 


programs in the area of stress analysis to obtain thermal stresses, 
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APPENDIX A 


In this section some solutions to different configurations will be 
presented as examples of possible problem types that may be solved with 
this program. Also the results of these examples are compared with 
analytic solutions of the same configurations. In order to avoid con- 
fusion the same material nickel (99.9%) and the same initial temperature, 
100°F, are used in each problem. 

The material properties of nickel are: 


e 


Cc 


556 1bm/ft3 


. 1065 BTU/1bm-°F 


il 


Xp = .882 ft2/hr 
Example 1. 

This is a problem in one dimensional heat conduction. It could also 
be considered a problem of heat conduction in a semi-infinite plate by 
assuming that lines parallel to the direction of heat flux are adiabatic 
boundaries. Because of this assumption we only need one row of elements 
across the plate which is one foot wide (see figure A-1). The top and 
bottom sides of the elements are adiabatic boundaries. Also the end 
boundary conditions are adiabatic on the left and are constant O°F on the 


right. The nodal points are numbered from left to right, starting with 


the top row. The elements are also numbered from left to right. The 


computer solution provides answers at distinct time intervals while most 





Figure A-1l One Dimensional Heat Conduction 
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analytic solutions are given in terms of the dimensionless time parameter 
Fourier number, Ngo. The Fourier number, Ngo, equals Bae where “fq 
is the thermal diffusivity of the material. The comparisons will be made 
at various Fourier numbers. The analytic solution was obtained from 


page 98 of reference (2]. 


Nodal NFO Carslaw and Finite 
Point Jaeger Element 
1 04 100 99.6 
5 04 99 99.0 
9 .04 96 95.9 
13 04 85 84.7 
17 . 04 54 54.5 
1 a! 95 94.8 
5 Pee 93 91.5 
9 wl 82 82.5 
13 1 63 63.6 
17 aul: 35 ia! 
1 4 48 48.1 
5 ae 45 45.8 
9 a 38 38.9 
13 4, 28 23), 3 
17 4, 15 14.9 

Example 2. 


This example has the same mesh configuration and boundary conditions 
as the first example. The initial temperature distribution is changed 
to the form T = 100 - 50x°F, The analytic solution was obtained from 


page 98 of reference (2]. 


Nodal Nero Carslaw and Finite 
Point Jaeger Element 
1 ~O1 95 94.0 
5 .O1 90 90.0 
9 .O1 80 79.8 
13 sO 1 70 68.0 
17 FOL 60 53.0 
1 04 99 89.0 
5 04 86 85.6 
9 04 77 76.9 
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13 . 04 62 G2. % 


17 04 36 37 
I I 80 80.0 
2, al 76 76.0 
9 sal 67 66.3 

13 I. 50 49,7 

OF, oll 28 27.0 

Example 3. 


This example is a solution for two-dimensional heat conduction. The 
configuration is a square plate with an initial temperature of 100°F and 
boundary conditions of O°F on all sides. Using symmetry the user need 


only examine one fourth of the plate, i.e. one corner (see figure A-2). 





Figure A-2 Two Dimensional Heat Conduction 


The nodal points are numbered from left to right and from the bottom row 
to the top. The elements are numbered in the same manner. In this case 
the problem consists of a plate with two sides at O°F and the other two 
insulated for an adiabatic boundary condition. The analytic solution is 
obtained by using the solution for one dimensional heat conduction in a 


semi-infinite plate and the method of multiplicative superposition. 


61 


The analytic solution was obtained from page 118 of reference (3]. 


Nodal NFO Finite 
Point Chapman Element 
8 O01 27.0 30.0 
15 O01 70.6 70.6 
22 O01 92.0 90.5 
29 .O1 97.0 9743 
36 .O1 100.0 98.8 
8 .05 5.75 6.29 
15 .05 2132 Wf ag 2 
22 .05 40.4 40.8 
29 205 54.0 DOG 
36 205 99.3 60.4 
8 aa): 1.96 2.23 
15 As | lee 3 Sel 
22 at 14.4 jay 72 
29 Fel 19.4 21.1 
39 a: 22.1 2353 

Example 4. 


In this problem the same mesh configuration as the previous example 
is used, The boundary conditions are changed so that the sides are adfa- 
batic, the top is at a constant 100°F. The initial temperature is 100°F, 
By taking an extremely large time increment, that is of a magnitude 
greater than the maximum number of significant figures available in the 
computer, say 1017, the first time step will yield the steady state 
solution to the problem. This is made possible because the time depend- 
ent contributions to the coefficient matrix and the forcing vector are 
multiplied by the inverse of the time increment which in this case is 
essentially zero (i97*7 5, The solution to this problem Js 4 linear 


temperature distribution and the results below were obt ajned by using a 


time increment of 1017 minutes, 


Nodal Finite 
Point Element. 
3 0.0 
9 20.0 
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KEY 
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LPRINT 
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MB AND 
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MTYPE 
NSEQ 
NUMCBC 
NUMEL 


NUMM AT 


APPENDIX B 


coefficient matrix 

dummy logical variable 

effective forcing vector 

ky in an element 

ky in an element 

kyy in an element 

convection boundary layer heat transfer coefficient 
heat generation in a material 

output print interval 

number of time steps 


k 


if 


1 - 4 nodal point numbers 


k 5 material number 


If 


temperature and flux boundary generation 
geometry Cype 

convection boundary generation 

temperature or flux boundary parameter 
spacial temperature distribution parameter 
print parameter | 
logical variable 

half band width 

logical variable 

material number 

number of time sequences 

number of convection boundary conditions 

number of elements 


number of materials 
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